{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%load_ext Cython"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%cython\n",
    "\n",
    "from libc.stdlib cimport malloc, free\n",
    "\n",
    "DEF LIMIT = 1024 * 31\n",
    "DEF PRIME = 1024 * 4\n",
    "DEF SIEVE = 1024 * 32\n",
    "\n",
    "cdef inline int imin(int a, int b) nogil:\n",
    "    return a if a < b else b\n",
    "\n",
    "cdef inline int memset(char *p, int n) nogil:\n",
    "    cdef:\n",
    "        short *q = <short *>p\n",
    "        int i, j = 0\n",
    "\n",
    "    for i in range((n + 1) >> 1):\n",
    "        j += q[i]\n",
    "        q[i] = 0x0100\n",
    "\n",
    "    return j >> 8\n",
    "\n",
    "cdef int naive_sieve(char *sieve, int *primes, int *offsets, int n) nogil:\n",
    "    cdef int i, j\n",
    "\n",
    "    memset(sieve, n)\n",
    "\n",
    "    for i in range(3, n, 2):\n",
    "        if sieve[i]:\n",
    "            j = i * i\n",
    "            while j < n:\n",
    "                sieve[j] = 0\n",
    "                j += i << 1\n",
    "\n",
    "            primes[0] = i\n",
    "            offsets[0] = j\n",
    "            primes += 1\n",
    "            offsets += 1\n",
    "\n",
    "    primes[0] = 0\n",
    "    offsets[0] = 0\n",
    "\n",
    "    return memset(sieve, n)\n",
    "\n",
    "cdef int segmented_sieve(char *sieve, int *primes, int *offsets, int k, int n) nogil:\n",
    "    cdef int i\n",
    "\n",
    "    while primes[0]:\n",
    "        i = offsets[0] - k\n",
    "        while i < n:\n",
    "            sieve[i] = 0\n",
    "            i += primes[0] << 1\n",
    "        offsets[0] = i + k\n",
    "\n",
    "        primes += 1\n",
    "        offsets += 1\n",
    "\n",
    "    return memset(sieve, n)\n",
    "\n",
    "cpdef int eratosthenes(int n) nogil:\n",
    "    cdef:\n",
    "        char *sieve\n",
    "        int *primes\n",
    "        int *offsets\n",
    "        int k, total\n",
    "\n",
    "    if n > LIMIT * LIMIT:\n",
    "        return -1\n",
    "\n",
    "    sieve = <char *>malloc(SIEVE)\n",
    "    primes = <int *>malloc(PRIME * sizeof(int))\n",
    "    offsets = <int *>malloc(PRIME * sizeof(int))\n",
    "\n",
    "    total = naive_sieve(sieve, primes, offsets, imin(n, LIMIT))\n",
    "\n",
    "    memset(sieve, SIEVE)\n",
    "    k = LIMIT\n",
    "    n -= LIMIT\n",
    "\n",
    "    while n > 0:\n",
    "        total += segmented_sieve(sieve, primes, offsets, k, imin(n, SIEVE))\n",
    "        k += SIEVE\n",
    "        n -= SIEVE\n",
    "\n",
    "    free(sieve)\n",
    "    free(primes)\n",
    "    free(offsets)\n",
    "\n",
    "    return total"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "primes below 10**1: 4\n",
      "primes below 10**2: 25\n",
      "primes below 10**3: 168\n",
      "primes below 10**4: 1229\n",
      "primes below 10**5: 9592\n",
      "primes below 10**6: 78498\n",
      "primes below 10**7: 664579\n",
      "primes below 10**8: 5761455\n",
      "primes below 10**9: 50847534\n"
     ]
    }
   ],
   "source": [
    "for i in range(1, 10):\n",
    "    print('primes below 10**%d: %d' % (i, eratosthenes(10**i)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## timeit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10000 loops, best of 3: 56.6 µs per loop\n",
      "1000 loops, best of 3: 574 µs per loop\n",
      "100 loops, best of 3: 6.07 ms per loop\n",
      "10 loops, best of 3: 68.4 ms per loop\n",
      "1 loop, best of 3: 863 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit eratosthenes(1024 * 31)\n",
    "%timeit eratosthenes(10**6)\n",
    "%timeit eratosthenes(10**7)\n",
    "%timeit eratosthenes(10**8)\n",
    "%timeit eratosthenes(10**9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
